1 Background

Considering the prevalence of obesity (1), its associated fracture risk [(2)] and that these patients are especially prone to develop loading associated musculoskeletal injuries [(5)], there is a need to develop ground reaction force (GRF) prediction equations that are accurate for overweight and obese subjects in order to precisely determine and monitor exercise-associated mechanical loading in these patients. Therefore, the purpose of this study was to develop prediction equations based on accelerometer data able to accurately predict GRF on a broad range of body masses, from normal weight to severely obese subjects, setting thereby, a base for objective prescription and monitoring of exercise mechanical loads.

2 Load data

Initially, 3 files needed to be loaded, each containing data from one of the accelerometer placements utilized in the study: ankle, back and hip.

A piece of a data frame used in all subsequent analysis can be seen below.

## # A tibble: 299 x 11
##       ID speed body_mass height   BMI sex     age pVGRF_N pVACC_g pRGRF_N
##    <dbl> <dbl>     <dbl>  <dbl> <dbl> <chr> <dbl>   <dbl>   <dbl>   <dbl>
##  1     3     2      85.6    171  29.3 M        50   1003.    2.21   1038.
##  2     3     3      85.6    171  29.3 M        50   1114.    2.84   1132.
##  3     3     4      85.6    171  29.3 M        50   1220.    3.64   1243.
##  4     3     5      85.6    171  29.3 M        50   1212.    3.78   1250.
##  5     3     6      85.6    171  29.3 M        50   1237.    3.76   1285.
##  6     4     2      61.1    162  23.3 F        45    670.    1.81    673.
##  7     4     3      61.1    162  23.3 F        45    694.    2.20    710.
##  8     4     4      61.1    162  23.3 F        45    743.    2.74    756.
##  9     4     5      61.1    162  23.3 F        45    803.    3.17    819.
## 10     4     6      61.1    162  23.3 F        45    872.    3.66    878.
## # … with 289 more rows, and 1 more variable: pRACC_g <dbl>

3 Ground reaction force X acceleration relationship

We first determined the relationship between pACC and pGRF obtained during the incremental walking speeds used in our experimental protocol. This was performed in order to verify the consistency of this relationship for increasing pACC values and the ability of the accelerometer to discriminate differences in pGRF between subjects in different BMI classes. A scatterplot with these relationships for all three accelerometer placements tested and for both the resultant and its vertical component is depicted in Figure 1. Generally, as expected, pACC recorded by accelerometers in all placements, showed a linear increase as a function of pGRF increases. Also, the recorded accelerations were shown to be able to provide a good discrimination between different BMI classes, as for the same registered pACC the pGRF tended to be consistently higher for subjects in higher BMI classes.

Code to generate these plots can be seen here.

4 Linear mixed models

Linear mixed models (LMM) were used to estimate peak resultant ground reaction force (pRGRF) and peak vertical ground reaction force (pVGRF). Distinct LMMs were developed with data from ankle, back and hip accelerometer placement (ankle, back and hip data frames) using the lme() function of the nlme package. Covariance structure used was an autoregressive process of order 1 (correlation = corAR1) and maximum likelihood method was used for estimating parameters (method = "ML"). Predictors tested on pRGRF models were body mass and peak resultant acceleration (pRACC), while on pVGRF were body mass and vertical acceleration (pVACC). All of them were tested as fixed effects and have shown to be significant predictors (e.g.: fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass). Both random intercept and slopes were tested, but only the random intercept inclusion has showed models improvement (random = ~ 1 | ID). Linear, quadratic and cubic polynomial simulations were also tested, whereas the last one did not contribute significantly to the models. Final models were chosen according to -2 log-likelihood statistics. Traditional coefficient of determination (R2) was represented by conditional R2 (6) computed with rsquared() function of the piecewiseSEM package, that estimates the variance explained by the whole model (7).

Code to build the LMMs, as well as their summary and R2, can be found below.

## Linear mixed-effects model fit by maximum likelihood
##  Data: ankle 
##        AIC      BIC    logLik
##   3188.582 3218.186 -1586.291
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02697863 86.14005
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8748119 
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       188.55057  62.95028 232  2.995230  0.0030
## pRACC_g            -2.04012  18.97302 232 -0.107527  0.9145
## I(pRACC_g^2)       -6.97317   2.21342 232 -3.150405  0.0018
## body_mass           6.56845   0.70675  62  9.293826  0.0000
## pRACC_g:body_mass   1.88108   0.15506 232 12.130984  0.0000
##  Correlation: 
##                   (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g           -0.700                      
## I(pRACC_g^2)       0.247 -0.700               
## body_mass         -0.898  0.426   0.130       
## pRACC_g:body_mass  0.668 -0.562  -0.170 -0.766
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.59736240 -0.63264008  0.03559703  0.69274354  2.95031219 
## 
## Number of Observations: 299
## Number of Groups: 64
##   Response   family     link method  Marginal Conditional
## 1  pRGRF_N gaussian identity   none 0.9297268   0.9297268
## Linear mixed-effects model fit by maximum likelihood
##  Data: back 
##        AIC      BIC    logLik
##   3185.997 3215.628 -1584.999
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:   0.0184402 79.53301
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8530994 
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       -698.7031 123.44832 233 -5.659884   0e+00
## pRACC_g           1047.5129 122.40109 233  8.558036   0e+00
## I(pRACC_g^2)      -345.2605  36.10073 233 -9.563809   0e+00
## body_mass            3.8294   1.06086  62  3.609729   6e-04
## pRACC_g:body_mass    6.0219   0.61893 233  9.729449   0e+00
##  Correlation: 
##                   (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g           -0.896                      
## I(pRACC_g^2)       0.648 -0.892               
## body_mass         -0.636  0.268   0.143       
## pRACC_g:body_mass  0.565 -0.281  -0.167 -0.922
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.0955967 -0.6001009 -0.0930397  0.5792458  3.1083008 
## 
## Number of Observations: 300
## Number of Groups: 64
##   Response   family     link method  Marginal Conditional
## 1  pRGRF_N gaussian identity   none 0.9361655   0.9361655
## Linear mixed-effects model fit by maximum likelihood
##  Data: hip 
##        AIC      BIC    logLik
##   3091.034 3120.637 -1537.517
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02080893  71.4172
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8662087 
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass 
##                       Value Std.Error  DF    t-value p-value
## (Intercept)       -300.9909  89.71869 232  -3.354829   9e-04
## pRACC_g            522.6850  71.61994 232   7.298037   0e+00
## I(pRACC_g^2)      -171.5606  16.91100 232 -10.144910   0e+00
## body_mass            3.9596   0.82234  62   4.814986   0e+00
## pRACC_g:body_mass    5.3671   0.41930 232  12.799992   0e+00
##  Correlation: 
##                   (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g           -0.878                      
## I(pRACC_g^2)       0.588 -0.854               
## body_mass         -0.766  0.422   0.030       
## pRACC_g:body_mass  0.678 -0.471  -0.038 -0.891
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.54273847 -0.59285086  0.05857753  0.70507501  2.51256437 
## 
## Number of Observations: 299
## Number of Groups: 64
##   Response   family     link method  Marginal Conditional
## 1  pRGRF_N gaussian identity   none 0.9488446   0.9488446
## Linear mixed-effects model fit by maximum likelihood
##  Data: ankle 
##        AIC      BIC    logLik
##   3190.571 3220.175 -1587.286
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02781982 91.52233
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8928218 
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       173.06487  77.90449 232  2.221500  0.0273
## pVACC_g            17.98949  39.07008 232  0.460442  0.6456
## I(pVACC_g^2)      -24.86971   6.71678 232 -3.702623  0.0003
## body_mass           5.35461   0.82437  62  6.495401  0.0000
## pVACC_g:body_mass   3.19203   0.27580 232 11.573831  0.0000
##  Correlation: 
##                   (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g           -0.767                      
## I(pVACC_g^2)       0.393 -0.785               
## body_mass         -0.856  0.414   0.081       
## pVACC_g:body_mass  0.680 -0.526  -0.092 -0.804
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.84130261 -0.57735799  0.02962816  0.55766690  2.61025073 
## 
## Number of Observations: 299
## Number of Groups: 64
##   Response   family     link method  Marginal Conditional
## 1  pVGRF_N gaussian identity   none 0.9177371   0.9177371
## Linear mixed-effects model fit by maximum likelihood
##  Data: back 
##        AIC      BIC    logLik
##   3234.998 3264.628 -1609.499
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02508531 102.5201
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.9079104 
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       -776.8934 145.57197 233 -5.336834       0
## pVACC_g           1042.9052 146.96000 233  7.096524       0
## I(pVACC_g^2)      -336.2115  44.82805 233 -7.500024       0
## body_mass            6.2132   1.22304  62  5.080120       0
## pVACC_g:body_mass    5.0805   0.72861 233  6.972936       0
##  Correlation: 
##                   (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g           -0.889                      
## I(pVACC_g^2)       0.652 -0.894               
## body_mass         -0.688  0.327   0.066       
## pVACC_g:body_mass  0.587 -0.341  -0.101 -0.891
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2287736 -0.6177471 -0.1216603  0.4585545  3.7585643 
## 
## Number of Observations: 300
## Number of Groups: 64
##   Response   family     link method  Marginal Conditional
## 1  pVGRF_N gaussian identity   none 0.8935621   0.8935621
## Linear mixed-effects model fit by maximum likelihood
##  Data: hip 
##        AIC      BIC    logLik
##   3153.414 3183.017 -1568.707
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept) Residual
## StdDev:  0.02202366 84.32584
## 
## Correlation Structure: AR(1)
##  Formula: ~1 | ID 
##  Parameter estimate(s):
##       Phi 
## 0.8869442 
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)       -435.7365 122.69794 232 -3.551294   5e-04
## pVACC_g            586.6627 112.66171 232  5.207295   0e+00
## I(pVACC_g^2)      -188.9689  28.55314 232 -6.618146   0e+00
## body_mass            5.8047   0.97485  62  5.954511   0e+00
## pVACC_g:body_mass    4.9544   0.54067 232  9.163472   0e+00
##  Correlation: 
##                   (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g           -0.910                      
## I(pVACC_g^2)       0.715 -0.907               
## body_mass         -0.759  0.471  -0.125       
## pVACC_g:body_mass  0.688 -0.524   0.135 -0.888
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.95040403 -0.53289689 -0.04360304  0.45305786  4.09646433 
## 
## Number of Observations: 299
## Number of Groups: 64
##   Response   family     link method  Marginal Conditional
## 1  pVGRF_N gaussian identity   none 0.9263472   0.9263473

In all models, conditional R2 values ranged from 0.89 to 0.95.

5 Validation analysis

Model validation was assessed by the leave-one-out cross-validation (LOOCV) method. This approach is recommended when a different sample is not available for cross-validation and it provides an insight on the model potential to predict outcomes in a new independent sample (8). For LOOCV each participant’s data was separated into a testing dataset (one at a time) with the remaining data being in the training dataset. New LMM, with the same outcomes and predictors as determined for the entire sample, were developed using the training dataset and then used to predict pGRF for the only participant in the testing dataset. This process was repeated for all participants (64 times). Data from testing dataset were used in the remaining statistical analysis.

LOOCV of the LMM was done with the cross_validate_mixed_model() function and is shown below.

5.1 Bland-Altman plots

Bland-Altman plots were used to examine the agreement between pGRF measured with force plates and those predicted through the regression equations. The difference of the actual and predicted pGRF was plotted against their mean. Bias was expressed as the mean of these differences and the limits of agreement were obtained using ±1.96 standard deviation of the mean between actual and predicted pGRF (9).

Bland-Altman plots can be seen below for all accelerometer placements. Panel A shows plots of pRGRF and panel B shows plots of pVGRF.

Code to generate these plots can be seen here.


One-sample t-tests were run to check whether bias from each accelerometer placement, in both pRGRF and pVGRF, was significantly different than 0. These tests were run using the t.test() function of the base R with the argument mu = 0. No significant differences were found.

## 
##  One Sample t-test
## 
## data:  LOOCV_ankle_res_LMM$diff
## t = 0.37558, df = 298, p-value = 0.7075
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -9.266825 13.638154
## sample estimates:
## mean of x 
##  2.185665
## 
##  One Sample t-test
## 
## data:  LOOCV_back_res_LMM$diff
## t = 0.44993, df = 299, p-value = 0.6531
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7.724802 12.303988
## sample estimates:
## mean of x 
##  2.289593
## 
##  One Sample t-test
## 
## data:  LOOCV_hip_res_LMM$diff
## t = 0.36236, df = 298, p-value = 0.7173
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7.774785 11.284133
## sample estimates:
## mean of x 
##  1.754674
## 
##  One Sample t-test
## 
## data:  LOOCV_ankle_vert_LMM$diff
## t = -0.77815, df = 298, p-value = 0.4371
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -17.18212   7.44448
## sample estimates:
## mean of x 
## -4.868821
## 
##  One Sample t-test
## 
## data:  LOOCV_back_vert_LMM$diff
## t = -0.10563, df = 299, p-value = 0.9159
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -12.19819  10.95537
## sample estimates:
##  mean of x 
## -0.6214066
## 
##  One Sample t-test
## 
## data:  LOOCV_hip_vert_LMM$diff
## t = 0.0087888, df = 298, p-value = 0.993
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -9.956701 10.046032
## sample estimates:
##  mean of x 
## 0.04466578

Also, linear regressions were applied to identify if there was any proportional bias, that is, if bias was related with the magnitude of the mean between measured and predicted pGRF (10). Linear regressions were run using the lm() function of the base R.

For lower back and hip placements, results showed a significant proportional bias (p < 0.05), however, with a low magnitude (highest R2 = 0.032). These results suggest that despite there is a trend for underestimation at increasingly higher pGRF values, the magnitude of this effect is neglectable. No proportional bias (p > 0.05) was detected for ankle derived equations.

## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_ankle_res_LMM)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -338.28  -54.43    3.59   57.89  285.63 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -33.03516   22.46494  -1.471    0.142
## mean          0.03011    0.01855   1.623    0.106
## 
## Residual standard error: 100.4 on 297 degrees of freedom
## Multiple R-squared:  0.00879,    Adjusted R-squared:  0.005453 
## F-statistic: 2.634 on 1 and 297 DF,  p-value: 0.1057
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_res_LMM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -267.995  -55.071   -2.505   52.450  294.243 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -47.77144   19.69289  -2.426  0.01587 * 
## mean          0.04282    0.01628   2.630  0.00899 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 87.28 on 298 degrees of freedom
## Multiple R-squared:  0.02268,    Adjusted R-squared:  0.0194 
## F-statistic: 6.915 on 1 and 298 DF,  p-value: 0.008991
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_res_LMM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -268.721  -53.909    5.937   51.546  235.762 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -40.38734   18.75117  -2.154   0.0321 *
## mean          0.03610    0.01552   2.325   0.0207 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 83.12 on 297 degrees of freedom
## Multiple R-squared:  0.01788,    Adjusted R-squared:  0.01457 
## F-statistic: 5.406 on 1 and 297 DF,  p-value: 0.02074
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_ankle_vert_LMM)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -418.05  -60.21   10.54   59.08  293.87 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -45.43603   24.52661  -1.853   0.0649 .
## mean          0.03523    0.02060   1.710   0.0883 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 107.8 on 297 degrees of freedom
## Multiple R-squared:  0.009752,   Adjusted R-squared:  0.006418 
## F-statistic: 2.925 on 1 and 297 DF,  p-value: 0.08827
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_vert_LMM)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -342.85  -64.14   -4.92   51.55  354.11 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -46.47386   23.16029  -2.007   0.0457 *
## mean          0.03998    0.01954   2.046   0.0416 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 101.4 on 298 degrees of freedom
## Multiple R-squared:  0.01386,    Adjusted R-squared:  0.01055 
## F-statistic: 4.187 on 1 and 298 DF,  p-value: 0.04162
## 
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_vert_LMM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -255.386  -53.520   -0.846   49.427  295.833 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -42.10744   19.93510  -2.112   0.0355 *
## mean          0.03683    0.01685   2.186   0.0296 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 87.33 on 297 degrees of freedom
## Multiple R-squared:  0.01583,    Adjusted R-squared:  0.01252 
## F-statistic: 4.778 on 1 and 297 DF,  p-value: 0.02961

5.2 Indices of accuracy

To evaluate models prediction accuracy, mean absolute error (MAE), mean absolute percent error (MAPE) and root mean square error (RMSE) were computed These indices were computed with the accuracy_indices() function.

Table below shows these indices for each accelerometer placement in both pRGRF and pVGRF.

Accelerometer Placement GRF MAE MAPE RMSE
Ankle pRGRF 77.0 7.0% 100.5
Ankle pVGRF 81.0 7.5% 108.1
Back pRGRF 66.7 5.9% 88.0
Back pVGRF 74.9 6.5% 101.7
Hip pRGRF 65.5 5.9% 83.6
Hip pVGRF 65.6 5.9% 87.7

All conditions showed a MAE from near 65 N to near 81 N, MAPE from near 6% to near 7% and RMSE from near 84 N to 84 N.

5.3 ANOVA

A series of repeated measures analysis of variance (ANOVA) were run to assess whether pGRF predicted from the regression equations were significantly different from those measured with FP. Walking speeds, accelerometer placements (ankle, back and hip), and the interaction effect (speed X accelerometer placements) were considered for analysis. These procedures were taken separately for resultant and its vertical component. If assumptions of sphericity were violated (p < 0.05), the conservative Greenhouse–Geisser correction factor would be applied to adjust the degrees of freedom. Post-hoc analyses would be conducted using pairwise comparisons with Holm’s test if a significant difference were observed among actual and predicted pGRF.

In order to run the repeated measures ANOVA, new data frames needed to be built, for both pRGRF and pVGRF, putting all the pGRF values in a single variable and grouping then by accelerometer placement or actual value in another variable. An example of code to tidy the data is shown below.

## Predicted pRGRF
ankle_res_pred <- LOOCV_ankle_res_LMM %>% 
  select(ID, speed, pRGRF_N_predicted) %>% 
  spread(key = speed, value = pRGRF_N_predicted) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = ankle
  )

back_res_pred <- LOOCV_back_res_LMM %>% 
  select(ID, speed, pRGRF_N_predicted) %>% 
  spread(key = speed, value = pRGRF_N_predicted) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = back
  )

hip_res_pred <- LOOCV_hip_res_LMM %>% 
  select(ID, speed, pRGRF_N_predicted) %>% 
  spread(key = speed, value = pRGRF_N_predicted) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = hip
  )

res_pred_df <- ankle_res_pred %>% 
  full_join(back_res_pred, by = c("ID", "speed")) %>% 
  full_join(hip_res_pred, by = c("ID", "speed")) %>% 
  na.omit()

## Actual pRGRF
ankle_res_actual <- LOOCV_ankle_res_LMM %>% 
  select(ID, speed, pRGRF_N) %>% 
  spread(key = speed, value = pRGRF_N) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_ankle
  )

back_res_actual <- LOOCV_back_res_LMM %>% 
  select(ID, speed, pRGRF_N) %>% 
  spread(key = speed, value = pRGRF_N) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_back
  )

hip_res_actual <- LOOCV_hip_res_LMM %>% 
  select(ID, speed, pRGRF_N) %>% 
  spread(key = speed, value = pRGRF_N) %>% 
  na.omit() %>% 
  gather(
    `2`, `3`, `4`, `5`, `6`,
    key = speed,
    value = actual_hip
  )

res_actual_df <- ankle_res_actual %>% 
  full_join(back_res_actual, by = c("ID", "speed")) %>% 
  full_join(hip_res_actual, by = c("ID", "speed")) %>% 
  na.omit()

res_actual_df <- res_actual_df %>% 
  mutate(actual = (actual_ankle + actual_back + actual_hip) / 3) %>% 
  select(ID, speed, actual) 

## Merge predicted and actual data frames
res_ANOVA_df <- res_pred_df %>% 
  full_join(res_actual_df, by = c("ID", "speed")) %>% 
  gather(
    ankle, back, hip, actual,
    key = "group",
    value = "pRGRF"
  )
res_ANOVA_df$ID    <- as.factor(res_ANOVA_df$ID)
res_ANOVA_df$speed <- as.factor(res_ANOVA_df$speed)
res_ANOVA_df$group <- as.factor(res_ANOVA_df$group)

Repeated measures ANOVAs were performed using the ezANOVA() function of the ez package. Summary statistics for pRGRF and pVGRF ANOVA can be seen below, respectively.

## $ANOVA
##        Effect DFn DFd          SSn        SSd            F             p
## 1 (Intercept)   1  45 1.121949e+09 55569779.5 9.085464e+02  1.758136e-31
## 2       speed   4 180 1.012027e+07   439209.6 1.036890e+03 4.502391e-123
## 3       group   3 135 1.815789e+03  1484985.0 5.502448e-02  9.829444e-01
## 4 speed:group  12 540 2.892448e+04   508964.5 2.557352e+00  2.676332e-03
##   p<.05          ges
## 1     * 9.508430e-01
## 2     * 1.485583e-01
## 3       3.130415e-05
## 4     * 4.984240e-04
## 
## $`Mauchly's Test for Sphericity`
##        Effect            W            p p<.05
## 2       speed 3.520032e-02 1.041884e-26     *
## 3       group 3.731595e-01 3.597677e-08     *
## 4 speed:group 3.580631e-06 1.091425e-64     *
## 
## $`Sphericity Corrections`
##        Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF]
## 2       speed 0.3944080 3.843045e-50         * 0.4063060 1.412612e-51
## 3       group 0.6892659 9.505875e-01           0.7229903 9.560913e-01
## 4 speed:group 0.2458901 5.887449e-02           0.2650287 5.417688e-02
##   p[HF]<.05
## 2         *
## 3          
## 4
## $ANOVA
##        Effect DFn DFd          SSn        SSd           F             p
## 1 (Intercept)   1  45 1.082599e+09 52600265.0 926.1727721  1.163820e-31
## 2       speed   4 180 7.572969e+06   364333.1 935.3627935 3.185492e-119
## 3       group   3 135 1.658433e+04  2164731.4   0.3447516  7.929979e-01
## 4 speed:group  12 540 2.257064e+04   474638.4   2.1399001  1.338656e-02
##   p<.05          ges
## 1     * 0.9511475606
## 2     * 0.1198692000
## 3       0.0002981690
## 4     * 0.0004057531
## 
## $`Mauchly's Test for Sphericity`
##        Effect            W            p p<.05
## 2       speed 4.088323e-02 2.251569e-25     *
## 3       group 6.550476e-01 2.393741e-03     *
## 4 speed:group 8.739178e-06 4.965719e-58     *
## 
## $`Sphericity Corrections`
##        Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF]
## 2       speed 0.4424565 3.134415e-54         * 0.4592394 3.449388e-56
## 3       group 0.7989465 7.473983e-01           0.8469999 7.594651e-01
## 4 speed:group 0.2757896 9.127102e-02           0.3002998 8.521643e-02
##   p[HF]<.05
## 2         *
## 3          
## 4

For either resultant and its vertical component, repeated measures ANOVA demonstrated that actual and predicted pGRF increased significantly (p < 0.001) along with increments in walking speed.


To assess for any significant differences between actual and predicted pGRF in each speed, five separate repeated measures ANOVA were done, one for each walking speed.

# For resultant peak ground reaction force
# 2 km/h
res_ANOVA_s2 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 2),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 3 km/h
res_ANOVA_s3 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 3),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 4 km/h
res_ANOVA_s4 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 4),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 5 km/h
res_ANOVA_s5 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 5),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 6 km/h
res_ANOVA_s6 <- ezANOVA(
  data     = res_ANOVA_df %>% filter(speed == 6),
  dv       = pRGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# For vertical peak ground reaction force
# 2 km/h
vert_ANOVA_s2 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 2),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 3 km/h
vert_ANOVA_s3 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 3),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 4 km/h
vert_ANOVA_s4 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 4),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 5 km/h
vert_ANOVA_s5 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 5),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

# 6 km/h
vert_ANOVA_s6 <- ezANOVA(
  data     = vert_ANOVA_df %>% filter(speed == 6),
  dv       = pVGRF,
  wid      = ID,
  within   = group,
  detailed = TRUE,
  type     = 3
)

No significant differences (p > 0.05) between actual and predicted pGRF in each speed were found.

These results can be seen in figure below.

6 Regressions comparison

Our equation was compared with a previously published reference equation (11), in which Neugebauer et al. used a similar approach for the prediction of pGRF. Comparison was performed using the equation to predict pVGRF from hip-worn accelerometers, as it was the only suitable with our data. The analysis was performed in three ways using: i) a subsample of participants with normal weight or overweight (BMI ≥ 18.5 and < 30 kg.m-2); ii) a subsample of obese participants (BMI ≥ 30 kg.m-2); iii) the whole sample.

New data frames needed to be built for each of the subsamples and the whole sample. Then, pVGRF was predicted using the reference equation. Code to build the data frames is shown below.

Bland-Altman plots were used to confront the agreement between pVGRF measured with FP and those predicted using both regression equations.

Code to generate these plots can be seen here.

A lower dispersion around the bias value, as well as a lower percentage of values falling off the LoA can be appreciated in our equation. Additionally, while the bias in our equation tended to be zero, in the reference equation bias was always > 0 (p < 0.05), showing a consistent underestimation of pVGRF.


Also, to assess the prediction accuracy, MAE, MAPE and RMSE were calculated using pVGRF values predicted from both equations.

The accuracy indices from our equation were all substantially lower compared to the reference equation indices, with an overall MAPE approximately 3 times smaller for our equation. For both equations, the MAPE was lower in the obese subsample than in the non-obese subsample

7 R session info

## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.5.2 (2018-12-20)
##  os       macOS Mojave 10.14.5        
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Maceio              
##  date     2019-07-18                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  abind          1.4-5    2016-07-21 [1] CRAN (R 3.5.0)
##  assertthat     0.2.0    2017-04-11 [1] CRAN (R 3.5.0)
##  backports      1.1.2    2017-12-13 [1] CRAN (R 3.5.0)
##  broman       * 0.68-2   2018-07-25 [1] CRAN (R 3.5.0)
##  broom          0.4.5    2018-07-03 [1] CRAN (R 3.5.0)
##  callr          2.0.4    2018-05-15 [1] CRAN (R 3.5.0)
##  car            3.0-2    2018-08-23 [1] CRAN (R 3.5.0)
##  carData        3.0-2    2018-09-30 [1] CRAN (R 3.5.0)
##  cellranger     1.1.0    2016-07-27 [1] CRAN (R 3.5.0)
##  cli            1.0.1    2018-09-25 [1] CRAN (R 3.5.0)
##  colorspace     1.3-2    2016-12-14 [1] CRAN (R 3.5.0)
##  cowplot      * 0.9.3    2018-07-15 [1] CRAN (R 3.5.0)
##  crayon         1.3.4    2017-09-16 [1] CRAN (R 3.5.0)
##  curl           3.2      2018-03-28 [1] CRAN (R 3.5.0)
##  data.table     1.12.0   2019-01-13 [1] CRAN (R 3.5.2)
##  desc           1.2.0    2018-05-01 [1] CRAN (R 3.5.0)
##  devtools       2.0.1    2018-10-26 [1] CRAN (R 3.5.0)
##  digest         0.6.19   2019-05-20 [1] CRAN (R 3.5.2)
##  dplyr        * 0.8.0.1  2019-02-15 [1] CRAN (R 3.5.2)
##  evaluate       0.10.1   2017-06-24 [1] CRAN (R 3.5.0)
##  ez           * 4.4-0    2016-11-02 [1] CRAN (R 3.5.0)
##  fansi          0.4.0    2018-10-05 [1] CRAN (R 3.5.0)
##  forcats      * 0.3.0    2018-02-19 [1] CRAN (R 3.5.0)
##  foreign        0.8-71   2018-07-20 [1] CRAN (R 3.5.2)
##  fs             1.2.6    2018-08-23 [1] CRAN (R 3.5.0)
##  ggplot2      * 3.1.0    2018-10-25 [1] CRAN (R 3.5.0)
##  glue           1.3.0    2018-07-17 [1] CRAN (R 3.5.0)
##  gtable         0.2.0    2016-02-26 [1] CRAN (R 3.5.0)
##  haven          2.0.0    2018-11-22 [1] CRAN (R 3.5.0)
##  here         * 0.1      2017-05-28 [1] CRAN (R 3.5.0)
##  highr          0.7      2018-06-09 [1] CRAN (R 3.5.0)
##  hms            0.4.2    2018-03-10 [1] CRAN (R 3.5.0)
##  htmltools      0.3.6    2017-04-28 [1] CRAN (R 3.5.0)
##  httr           1.3.1    2017-08-20 [1] CRAN (R 3.5.0)
##  jsonlite       1.5      2017-06-01 [1] CRAN (R 3.5.0)
##  knitr        * 1.21     2018-12-10 [1] CRAN (R 3.5.0)
##  labeling       0.3      2014-08-23 [1] CRAN (R 3.5.0)
##  lattice        0.20-38  2018-11-04 [1] CRAN (R 3.5.2)
##  lazyeval       0.2.1    2017-10-29 [1] CRAN (R 3.5.0)
##  lme4           1.1-19   2018-11-10 [1] CRAN (R 3.5.0)
##  lubridate      1.7.4    2018-04-11 [1] CRAN (R 3.5.0)
##  magrittr       1.5      2014-11-22 [1] CRAN (R 3.5.0)
##  MASS           7.3-51.1 2018-11-01 [1] CRAN (R 3.5.2)
##  Matrix         1.2-15   2018-11-01 [1] CRAN (R 3.5.2)
##  memoise        1.1.0    2017-04-21 [1] CRAN (R 3.5.0)
##  mgcv           1.8-26   2018-11-21 [1] CRAN (R 3.5.2)
##  minqa          1.2.4    2014-10-09 [1] CRAN (R 3.5.0)
##  mnormt         1.5-5    2016-10-15 [1] CRAN (R 3.5.0)
##  modelr         0.1.2    2018-05-11 [1] CRAN (R 3.5.0)
##  munsell        0.5.0    2018-06-12 [1] CRAN (R 3.5.0)
##  nlme         * 3.1-137  2018-04-07 [1] CRAN (R 3.5.2)
##  nloptr         1.2.1    2018-10-03 [1] CRAN (R 3.5.0)
##  openxlsx       4.1.0.1  2019-05-28 [1] CRAN (R 3.5.2)
##  pbkrtest       0.4-7    2017-03-15 [1] CRAN (R 3.5.0)
##  piecewiseSEM * 2.0.2    2018-07-24 [1] CRAN (R 3.5.0)
##  pillar         1.3.1    2018-12-15 [1] CRAN (R 3.5.0)
##  pkgbuild       1.0.2    2018-10-16 [1] CRAN (R 3.5.0)
##  pkgconfig      2.0.2    2018-08-16 [1] CRAN (R 3.5.0)
##  pkgload        1.0.2    2018-10-29 [1] CRAN (R 3.5.0)
##  plyr           1.8.4    2016-06-08 [1] CRAN (R 3.5.0)
##  prettyunits    1.0.2    2015-07-13 [1] CRAN (R 3.5.0)
##  processx       3.1.0    2018-05-15 [1] CRAN (R 3.5.0)
##  psych          1.8.4    2018-05-06 [1] CRAN (R 3.5.0)
##  purrr        * 0.3.0    2019-01-27 [1] CRAN (R 3.5.2)
##  R6             2.4.0    2019-02-14 [1] CRAN (R 3.5.2)
##  Rcpp           1.0.0    2018-11-07 [1] CRAN (R 3.5.0)
##  readr        * 1.3.0    2018-12-11 [1] CRAN (R 3.5.0)
##  readxl         1.1.0    2018-04-20 [1] CRAN (R 3.5.0)
##  remotes        2.0.2    2018-10-30 [1] CRAN (R 3.5.0)
##  reshape2       1.4.3    2017-12-11 [1] CRAN (R 3.5.0)
##  rio            0.5.16   2018-11-26 [1] CRAN (R 3.5.0)
##  rlang          0.3.1    2019-01-08 [1] CRAN (R 3.5.2)
##  rmarkdown      1.10     2018-06-11 [1] CRAN (R 3.5.0)
##  rprojroot      1.3-2    2018-01-03 [1] CRAN (R 3.5.0)
##  rstudioapi     0.7      2017-09-07 [1] CRAN (R 3.5.0)
##  rvest          0.3.2    2016-06-17 [1] CRAN (R 3.5.0)
##  scales         1.0.0    2018-08-09 [1] CRAN (R 3.5.0)
##  sessioninfo    1.1.1    2018-11-05 [1] CRAN (R 3.5.0)
##  stringi        1.3.1    2019-02-13 [1] CRAN (R 3.5.2)
##  stringr      * 1.4.0    2019-02-10 [1] CRAN (R 3.5.2)
##  testthat       2.0.0    2017-12-13 [1] CRAN (R 3.5.0)
##  tibble       * 2.0.1    2019-01-12 [1] CRAN (R 3.5.2)
##  tidyr        * 0.8.2    2018-10-28 [1] CRAN (R 3.5.0)
##  tidyselect     0.2.5    2018-10-11 [1] CRAN (R 3.5.0)
##  tidyverse    * 1.2.1    2017-11-14 [1] CRAN (R 3.5.0)
##  usethis        1.4.0    2018-08-14 [1] CRAN (R 3.5.0)
##  utf8           1.1.4    2018-05-24 [1] CRAN (R 3.5.0)
##  withr          2.1.2    2018-03-15 [1] CRAN (R 3.5.0)
##  xfun           0.4      2018-10-23 [1] CRAN (R 3.5.0)
##  xml2           1.2.0    2018-01-24 [1] CRAN (R 3.5.0)
##  yaml           2.1.19   2018-05-01 [1] CRAN (R 3.5.0)
##  zip            1.0.0    2017-04-25 [1] CRAN (R 3.5.0)
## 
## [1] /Library/Frameworks/R.framework/Versions/3.5/Resources/library

References

1. Guthold R, Stevens GA, Riley LM, Bull FC. Worldwide trends in insufficient physical activity from 2001 to 2016: a pooled analysis of 358 population-based surveys with 19 million participants. Lancet Glob Heal. 2018;6(10):e1077–86.

2. Li X, Gong X, Jiang W. Abdominal obesity and risk of hip fracture: a meta-analysis of prospective studies. Osteoporos Int. 2017;28(10):2747–57.

3. Hind K, Rudman H, Pearce M, Treadgold L, Birrell F. Obesity, bone density for weight and prevalent vertebral fracture at age 62 years: the Newcastle Thousand Families Study. J Clin Densitom. 2018;21(4):606.

4. Kaze AD, Rosen HN, Paik JM. A meta-analysis of the association between body mass index and risk of vertebral fracture. Osteoporos Int. 2018;29(1):31–9.

5. Viester L, Verhagen EA, Hengel KM, Koppes LL, Van Der Beek AJ, Bongers PM. The relation between body mass index and musculoskeletal symptoms in the working population. BMC Musculoskelet Disord. 2013;14(238).

6. Nakagawa S, Schielzeth H, O’Hara RB. A general and simple method for obtainingR2from generalized linear mixed-effects models. Methods in Ecology and Evolution. 2013;4(2):133–42.

7. Jaeger BC, Edwards LJ, Das K, Sen PK. An r2 statistic for fixed effects in the generalized linear mixed model. Journal of Applied Statistics. 2016;44(6):1086–105.

8. Staudenmayer J, Zhu W, Catellier DJ. Statistical considerations in the analysis of accelerometry-based activity monitor data. Medicine & Science in Sports & Exercise. 2012;44(1 Suppl 1):S61–7.

9. Bland JM, Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet. 1986;1(8476):301–10.

10. Giavarina D. Understanding bland altman analysis. Biochemia Medica. 2015;25(2):141–51.

11. Neugebauer JM, Collins KH, Hawkins DA. Ground reaction force estimates from ActiGraph GT3X+ hip accelerations. PLoS One. 2014/06/11. 2014;9(6):e99023.